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Abstract 

Using global interpolation functions (GIFs), 
boundary element solutions are obtained for two- 
and three-dimensional viscous flows. The so- 
lution is obtained in the form of a boundary 
integral plus a series of global basis functions. 
The unknown coefficients of the GIFs are deter- 
mined to ensure the satisfaction of the govern- 
ing equations at selected collocation points. The 
values of the coefficients involved in the bound- 
ary integral equations are determined by enforc- 
ing the boundary conditions. Both primitive- 
variable and vorticity-velocitv formulations are 
examined. 


[1993]: Lafe [1993]: Lafe £- Cheng [1994]). the 
method is being proposed for certain classes of 
nonlinear problems. 

Using the Dual Reciprocity approach, a given 
problem is typically decomposed into two parts 
- the linear and nonlinear portions. The solution 
to the linear portion is represented by a bound- 
ary integral whose kernel consists of the funda- 
mental solution to the linear governing equation. 
The nonlinear part is represented by either 1) lo- 
cal basis functions ( Brebbia et a/., [1991]): or 
2) global interpolation functions (GIFS) (Lafe 
[1993]). In either case, the boundary integral 
expressions and interpolation functions contain 
coefficients whose values are to be determined 
by enforcing the boundary conditions. When the 
“direct BEM” approach is followed the unknown 
coefficients are in essence the unknown physical 
variables (velocity components, pressure, tem- 
perature) of the problem. On the other hand, 
using the “indirect BEM” approach, the un- 
knowns are the weights/strengths of the bound- 
ary sources/dipoles and the local/global interpo- 
lating functions. The computational intensity of 
the indirect approach is much less than for the 
direct. 

In this paper, we report the formulation and 
development of GIF-based indirect BEM codes 
for two- and three-dimensional incompressible 
viscous flows. 


Governing Equations 

Introduction 2D Primitive Variable Formulation 


The boundary element method (BEM) has tra- 
ditionally been applied to problems governed by 
linear differential equations. At the core of the 
basic BEM computational process is the fun- 
damental solution (also referred to as the free- 
space Green’s function) defined as the impulse 
response of the governing equation to a unit ac- 
tion. This fundamental solution is either too dif- 
ficult or impossible to derive for practical non- 
linear problems. Recently, with the introduc- 
tion of the so-called Dual Reciprocity techniques 
(see e.g., Nardini & Brebbia [1982]; Brebbia et 
al., [1991]; Partridge et al., [1992]; Cheng et al., 
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where (u„,t>.) are the velocity components, 
and p. is the pressure. Let x = x./L; y = y./L; 



u = u./v m : v = v./v.; p = p./(pr„ 2 ). With 
these the governing equations become: 
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where R e = pv^Lfp. 

By taking the divergence of the momentum 
equation and *using the continuity equation, we 
obtain a PDE for the pressure in the form: 


We distribute sources (plus vortices * for the 
velocities) of unknown strength ^ on the flow 
boundary T. The combination of sources and 
vortices for the velocity components is important 
because of the need to automatically conserve 
mass. The solution for u % v,p is written as: 

«( x ) = J {wi( x, )ffii( x 4 x/ ) + w 2 (x')5 12 (x.x')} dT 

+ 5Z {^l*$*ll( x ) + $ 2 *$*12( x )} (8) 

* 

”( x ) = ^.{wi(x , )p2i( x ,x , ) + w 2 (x , )g 2 2( x .x / )} dr 

+ X] {/^l*$*2l( x ) + /?2*$*22( x )} (9) 

k 

p{x) = [ h*(x')033(x.x') dr + ]T#3*$*33 (10) 

JV k 


^ = 2 (l 5 -sb) < 7 > 

BEM-GIF Formulation The solution for a 
flow variable £ is decomposed into two parts (Co 
and Ci), such that 


C = Co + Cl 

and Co is the solution to the convection-free prob- 
lem 

v 2 Co = 0 

Therefore, the correction Ci is the addition re- 
quired such that the full governing equation is 
satisfied. The part Co can be modeled to a high 
degree of precision by boundary integral equa- 
tions. The correction Ci is represented by a series 
formed by globed basis functions. 

The fundamental solution for convection-free 
two-dimensional must satisfy 


in which x = (x,y), gij are associated with 
the fundamental solutions of the convection-free 
problem, and $*,y are global interpolation bases. 
For mass to be automatically conserved the free- 
space functions and GIFs must satisfy the fol- 
lowing conditions: 
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522 = g = In r, we have 
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V 2 y(x,x') = 2x<5(x,x') 

where 6 is the Dirac delta function applied at 
point x = (x,y), and felt at the point x = 
(x', y > ). The closed form solution for g is eas- 
ily shown to be (Jaswon & Symm [1977]) 

5 = In r 


Note that when gu is a source, the associated 
functions required to ensure mass conservation 
9ij ( l / j) turn out to be vortices. 

Similarly, using hyperbolic cosine global bases, 
with $*ii = $*22 = cosh(m*x)cosh(n*y) 
(m*,n* = 1,2.3,---) we have 

Tlfe 

$*12 = sinh(m*x)sinh(n*y) 

m* 

$*2i = — — - sinh(m*x)sinh(n*y) 


where r = |x — x'|. 
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The solutions for other auxiliary bases y^ij (i ^ 
j) for given are shown below: 

• Polynomial GIF 

*kn = $*22 = x m y n 

*kl2 = * m “V +1 

n + 1 

^21 = 2_x m+, y w " 1 

m + 1 

• Trigonometric GIF 


BEM-GIF Equations 

We distribute sources of strength (for ) and 
u >2 (for Q on the boundary T. The complete 
solution can be written in the form: 

$(x) = [u> 1 (x')dx' + ^/3 u .$*(x)(17) 

Jr k 

C(x) = f u)i{x')dx' + '*T 02k* k{*){ 18) 

Jr V 


= $*22 
*k 12 

$*21 

Radial GIF 
$*11 = $*22 = 

$*12 = 


cos(mx) cos(ny) 

— sin(mx)sin(ny) 
m 

— sin(mx)sin(ny) 
n 


— m = 2(n + 1) 
m 


The coefficients (i = 1,2) are determined by 
enforcing equations (15) and (16) at select collo- 
cation points, while w; (t = 1,2) are determined 
by enforcing the boundary conditions. 

3D Primitive Variable Formulation 


f he governing equations are 
du dv dw 

$*2 1 = -r"{(*-x n )(y-y n ) + U + ^ + ** ~ ° 
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If rid terms axe selected in the GIF series, then 
the coefficients /?;* (i = 1,2,3; k = 1,2 , •♦•nj) 
are determined by enforcing the momentum 
equations at nj collocation points. The strengths 
u (i = 1,2,3) of the sources /vortices are deter- 
mined by enforcing the boundary conditions at 
selected 74 boundary nodes. The details of the 
numerical implementation is outlined later. 
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2D Streamfunction-Vorticity Formulation 



In two-dimensional flow, the streamfunction $ 
and vorticity ( are defined by 

d$ 

U ~ dy 

dy 

v dx 

_ du dv 

dy dx 

When these are used in the governing equations, 
the continuity equation is automatically satis- 
fied, and the momentum equations become: 

V 2 * = -C (15) 

-*(£ + $) < 16 > 


Taking the divergence of the momentum equa- 
tion and enforcing mass conservation, the perti- 
nent PDE for the pressure is: 

■ -■{(£)(£)♦ (£)(£)*(£)(£)) 
“{(sewatsMtKS)) 

(23 


BEM-GIF Equations 

The fundamental solution for convection-free 
three-dimensional must satisfy 

V 2 <j(x,x') = 4tt£(x,x') 
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where 6 is the Dirac delta function applied at 
point x = (x,y.z), and felt at the point x = 
(x\ y \ z 1 ). The closed form solution for g is easily 
shown to be (J as won <L* Symm [1977]) 


5 = 


1 

r 


where r = |x — x'|. 

The full solution can be written in the form 

u(x) = J^ui{x')g u {x,x') dT 

+ J wa(x')p«(x,x') dr 

+ f u 3 (x')g 13 (x,x’) dT 
■ Jr 

+ ^20ik^k it (x) 

k 

+ £jS 2fc^fcl2(x) 

k 

+ 5Z /?3Jt^Jfcl3( x ) (24) 
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Sources/Dipoles By selecting yn = g?: = 
533 = 544 = g = 1/r and enforcing continuity, 
the functions g tJ (i ^ j) are given by: 


where 


and 
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u(x) = ^w 1 (x')52i(x.x') dr 
+ ^w 2 (x , ) 522 (x,x') dr 
+ ^w 3 (x') 523 (x,x') dr 

+ XI / ? lJt^Jfe2l(x) 

k 

+ 5Z^2fc^*22(x) 

k 

+ 5I# 3 *** 23 ( X ) ( 25 ) 

it 

«’(x) = ^u; 1 (x / )53i(x,x / ) dr 
+ jU(x')532(x.x') 

+ ^W3(x / )y 33 (x,x / ) dr 

+ 5Z ^lt^*3l(x) + /?2*$*32(x) 

Jt 
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ry-- = ( 5 -y , ) 2 + (---^) 2 

GIF Bases If $*11 = $*22 = $*33 = $*44 = 
cosh(mjtx) cosh(njty) cosh(/jtz) then for mass to 
be automatically conserved, the auxiliary func- 


tions $kij 

( i ^ j ) are given by: 

$*12 = 

n k 

. sinh(mjtx) sinh(njty) cosh(/* z) 

2m k 

4**13 = 

l k 

2m* sm M™Jt*)cosh(7i*y)sinh(/*-) 
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2n* 
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Jl k 
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Similar auxiliary functions can be derived for 
other selections of GIF such as 


p(x) = ^w 4 (x / )544(x,x / ) dr 

+ yi /^4* $*44 
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2. Trigonometric cos(mx) cos(ny) cos(lz); and 

3. Radial 1 + r 2 ( m+1 ) functions. 



3D Velocity- Vorticity Formulation BEM/GIF Equations 


The vorticity vector = (Cx-Cy-G) > s defined by We write the solution in the form: 


( = Vxu 

where Vx is the curl operator. Using the above 
in conjunction with the continuity equation V • 
u = 0 the following Poisson equations can be 


derived for the velocity components: 


V 2 u =dCy /dz-dCjdy 

(28) 

V 2 v = d<; z /dx - d( x /dz 

(29) 

V 2 tn = dCx/dy - d^y/dx 

(30) 

where x = (x,y, z ) is the spatial coordinate. 

The three-component vorticity 
equation is ( Fletcher [1991]) 

transport 

V-(uC)-(C-V)u-i-V 2 C = 

0 (31) 

For three-dimensional flows, there are six 
equations for it, t>,u?,Cr>Cy»G- These equations 
are expressible in the general elliptical form: 

V 2 0 = f 

(32) 


where 0 = #i = u, 9 2 = v, 

$3 — u. * 64 — Cx i ^5 — Cy? ^6 ~- G> and 
f = with: 


h I 


r ^Cv ftp ^ 

bz by 
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u 


Re (£ll ~ £12) 

/ 5 


Re (£21 — 622) 

. U . 


. Re (£31 — £32) j 


0i(x) = / w,(x')y(x.x')dx' + 5^ Jijt# fc (x) 

Jr k 

i = 1 , 2 , •• - 6 ( 34 ) 

where g is the free-space Green's function for the 
Laplace’s equation (i.e., convection-free flow). 

The coefficients are determined by enforc- 
ing the governing equations at select collocation 
points, while the strengths u;,- of the fictitious 
boundary sources are determined by enforcing 
the boundary conditions. 

Numerical Implementation 

Regardless of which formulation is used, the 
above BEM-GIF schemes share a common pat- 
tern: each solution consists of a sum of a bound- 
ary integral and a finite series. The numerical 
implementation can also be generalized for all 
the formulations. 

Consider the flow variable 6. We subdivide 
the boundary into elements. Let Nj(x) (j = 
1, 2, • • • n(,) represent the basis functions describ- 
ing the distribution of ui on T. 

Source/Vortex Strength Determination 

By selecting each of the nj, boundary points as 
successive origins of integration, the pertinent in- 
tegral equations can be assembled into the sys- 
tem: 
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71 fjj 

y! ) ’ a ikj LJ kj — i — 1, 2, • • • , nj, (35) 
j= 1 1 

where n w is the number of source/ vortex combi- 
nations used to model the flow variable, and 


&ikj 


bi 


Jr } Wj( x ')ffijy(x\xi) dx! x, e r e 
< (36) 

Ur,^(*0%^(x%x,)dx' x iG r, 

( n*i) - Eji, E"=1 to*ikj X, € T e 
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l X,- e r. 
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where Tj is the boundary segment on which the 
variable 6 is prescribed while T, is the portion 
on which d6/dn = V0 • n is prescribed. 

Therefore, we have equations to deter- 
mine Ujk (j = 1,2, •••nj k = 1,2, •••n u/ ). Sym- 
bolically equation (35) can be written in the al- 
ternative form: 

[A]{W} = (B} (38) 

which can be inverted to give: 

{W} = M]- 1 {6} (39) 

GIF Coefficients Determination 

The algebraic equations for determining the GIF 
coefficients 0 are nonlinear, on account of the 
convective terms. Therefore, the solution pro- 
cess has to be iterative. We start by assuming a 
convection-free flow, so that all the 0 values are 
set to zero initially. We then compute the error 
s in the transport equation for 0 at each of the 
714 collocation points. The error A 0 is assumed 
to be related to e through 

A0 = Xs 

where A is a relaxation parameter. The error A0 
is a function of the errors in the GIF coefficients. 
That is: 

n d n us 

^ = ( 4 °) 

k=lj=l 

Hence, by computing A0 for all nj points and all 
flow variables we have a system of equations: 

[*]{A/?} = {Ae} (41) 

which when inverted gives: 

{A/?} = f*]- 1 {\s} (42) 

The iterative steps are: 

1. Start with a trial (3k (k = 1,2, • • • n d ) for all 
flow variables. We have found a zero initial 
value to be the most convenient. 

2. Obtain W using equation (39). 


3. Use discretized forms of the appropriate in- 
tegral equations to compute 6. V0 at all n.j 
points. 

4. Compute the error s in the transport equa- 
tion at each of the collocation points. 

5. Use equation (42) to compute new values of 

P- 

6. Go back to Step 2 if convergence condition 
is still unsatisfied. 

Note that the matrix inversions in equations (39) 
and (42) need only be performed once, for fixed 
boundary problems. The vectors W and j3 are 
the quantities whose values change during the it- 
erative process. Once convergence is reached, the 
discretized integral equations can be used rou- 
tinely to obtain 0 = (u,v, tr,p, $,(,*. CyiG) or 
the gradient at any point (x) of interest. 


BEM Codes &: Preliminary Tests 

Four boundary element codes were developed. 

These include: 

1. PRIM-2D This is a two-dimensional gen- 
eral purpose boundary element solver based 
on the use of primitive variables. The code 
can be applied to any geometry. Bench- 
mark tests used global interpolation func- 
tions drawn from families of trigonometric, 
hyperbolic, radial, polynomial, wavelet, and 
Chebychev functions. 

2. VORT-2D This a two-dimensional code 
based on a streamfunction-vorticity formu- 
lation. The GIF bases enumerated in (1) 
were tested. 

3. PRIM-3D This a full three-dimensional 
boundary element code using primitive vari- 
ables. Four momentum equations are solved 
for the three velocity components (u,u, u/) 
and the pressure (p). 

4. VORT-3D A full three-dimensional code 
based on the vorticity- velocity formulation. 
Six transport equations are solved for the 



three vorticitv components ( Cx • Cy • G ) and 
the three velocity components (u.r.u ). 

Typical performance test results using PRIM- 
2D for the lid-driven cavity problem (Fig. 1) 
are shown in Figs. 2-3. The number of bound- 
ary elements «(, = 64 ( i.e 8 per side); and the 
number of collocation points nj = 64. The test 
results shown were obtained using the following 
relaxation parameters: 

A u = 0.2 
A p = 0.02 
A w = 1.0 

Preliminary test runs on three-dimensional de- 
veloping flows in straight and curved square 
ducts (Fig. 4) were also carried out. We set 
a = b = 1; R = 10, and a = 7 t/ 2. We assumed 
a plug flow (tz = 1; v = w — 0) at the inlet. A 
total of 72 rectangular boundary elements and 
343 internal nodes were used. Other salient pa- 
rameters include: 


L. and other parameters that may influence 
the flow problem. The convergence rate is 
highly sensitive to the choice of A. 

2. The optimal location of the collocation 
points Xfc (k = 1.2. •••nj) so that the en- 
suing solution process will capture the un- 
derlying physics of the problem. It is clear 
that nj can be much smaller than the num- 
ber of computational nodes required in the 
domain methods ( e.g ., finite elements and 
finite differences). However, since we are at 
liberty to select these points, their optimal 
location is crucial, in order to exploit the 
inherent advantages of the BEM-GIF for- 
mulation. 
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A u — A r — A u , — A w — 0.001 
A p = 1.0 
R c = 1000 

The velocity convergence characteristics of 
VORT-3D for the curved square duct is shown 
in Fig. 5. 

This is an ongoing research effort. We plan 
to report the comprehensive performance char- 
acteristics of all four BEM modules in an upcom- 
ing publication. 


Conclusions 


The paper presents the general foundation for 
obtaining boundary integral solutions to incom- 
pressible viscous flows. The GIF- based bound- 
ary integral equations ensure the automatic con- 
servation of mass. Critical factors which affect 
the convergence rates and the quality of the so- 
lutions include: 

1. The choice of the relaxation parameter A as 
a function of the Reynold’s number R e , the 
number of boundary elements nt, the num- 
ber of collocation points nj, the domain size 
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Figure 1: Lid-driven Cavity Problem 
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Figure 2: Velocity Convergence Rate for Cavity Problem 
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Figure 4: Curved Square Duct Problem 
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Figure 5: Typical Velocity Convergence Rate for Curved Square Duct 



